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Many physical processes, including the intensity fluctuations of a 
chaotic laser, the detection of single photons, and the Brownian mo¬ 
tion of a microscopic particle in a fluid are unpredictable, at least 
on long timescales. This unpredictability can be due to a variety of 
physical mechanisms, but it is quantified by an entropy rate. This 
rate, which describes how quickly a system produces new and random 
information, is fundamentally important in statistical mechanics and 
practically important for random number generation. We experimen¬ 
tally study entropy generation and the emergence of deterministic 
chaotic dynamics from discrete noise in a system that applies feed¬ 
back to a weak optical signal at the single-photon level. We show 
that the dynamics transition from shot noise to chaos as the photon 
rate increases, and that the entropy rate can reflect either the deter¬ 
ministic or noisy aspects of the system depending on the sampling 
rate and resolution. 

entropy | chaos | nonlinear dynamics | photon counting | statistical physics 

C ontinuous variables and dynamical equations are often 
used to model systems whose time evolution is comprised 
of discrete events occurring at random times. Examples in¬ 
clude the flow of ions across cell membranes , the dynamics 
of large populations of neurons [^, the birth and death of 
individuals in a population [^, traffic flow on roads 1^, the 
trading of securities in financial markets [^|^, infection and 
transmission of disease [^, and the emission and detection of 
photons . We can identify two sources of unpredictability in 
these systems: the noise associated with the underlying ran¬ 
dom occurrences that comprise these signals, which is often 
described by a Poisson process, and the macroscopic dynam¬ 
ics of the system, which may be chaotic. When both effects are 
present, the macroscopic dynamics can alter the statistics of 
the noise, and the small-scale noise can in turn feed the large- 
scale dynamics. This can lead to subtle and non-trivial effects 
including stochastic resonance and coherence resonance [ ^11| . 
Dynamical unpredictability and complexity are quantineoby 
Lyapunov exponents and dimensionality, while shot noise is 
characterized by statistical metrics like average rate, variance, 
and signal-to-noise ratio. Characterizing the unpredictability 
of a system with both large-scale dynamics and small-scale 
shot noise remains an important challenge in many disciplines 
including statistical mechanics and information security. 

Many cryp tographic applications, including public key en¬ 
cryption employ random numbers. Because the unpre¬ 
dictability of these numbers is essential, physical p rocesse s 
are sometimes used as a source of random numbers 
Physical r and om number ge ner ators are usually tested using 
the NIST and Diehard test suites which assess their 
ability to produce bits that are free of bias and correlation. 
These tests are an excellent assessment of the performance of a 
physical random number generator in practical situations, but 
leave an important and fundamental problem unaddressed. 
Deterministic post-processing procedures, such as hash func¬ 
tions are often employed to to remove bias and correlation. 
Because these procedures are algorithmic and reproducible. 


they cannot in principle increase the entropy rate of a bit 
stream. Thus, the reliability of a physical random number 
generator depends on an accurate assessment of the entropy 
rate of physical process that generated the numbers . It re¬ 
mains difficult to assess the unpredictability of a system based 
on physical principles. 

Evaluation of entropy rates from an information-theoretic 
per spectiv e is also centrally important in statistical mechan¬ 
ics [3b}f37] . One might expect that the unpredictability of a 
system with both small-scale shot noise and large-scale chaotic 
dynamics would depend on the scale at which it is observed. In 
many systems, the dependence of the entropy rate on the reso¬ 
lution, £, and the samplin g int erval, r, can reflect the physical 
origin of unpredictability [38}{41| . This dependence has been 
been studied experimentall^nBrownian motion, RC circuits, 
and Rayleigh-Benard convection [3 ^3^[3^[42 4^ . 

Here, we present an experimental exploration and numer¬ 
ical model of entropy production in a photon-counting opto¬ 
electronic feedback oscillator. Optoelectronic feedback loops 
which employ analog detectors and macroscopic optical sig¬ 
nals produce rich dyna mics wh ose timescales and dimension¬ 
ality are highly tunable E 48 . Our system applies optoelec¬ 
tronic feedback to a weak optical signal which is measured by 
a photon-counting detector. The dynamic range of this system 
(eight orders of magnitude in time scale and a factor of 256 in 
photon rate) allows us to directly observe the transition from 
shot-noise dominated behavior to a low-dimensional chaotic 
attractor with increasing optical power - a transition which 


Significance 

The unpredictability of physical systems can depend on the scale 
at which they are observed. For example, single photons inci¬ 
dent on a detector arrive at random times, but slow intensity 
variations can be observed by counting many photons over large 
time windows. We describe an experiment in which we modulate 
a weak optical signal using feedback from a single-photon de¬ 
tector. We quantitatively demonstrate a transition from single¬ 
photon shot noise to deterministic chaos. Furthermore, we show 
that measurements of the entropy rate of a system with small 
scale noise and large scale deterministic fluctuations can resolve 
both behaviors. We describe how quantifying entropy produc¬ 
tion can be used to evaluate physical random number generators. 
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to our knowledge has never been observed experimentally. We 
show that the entropy rate can reflect either the deterministic 
or stochastic aspects of the system, depending on the sampling 
rate and measurement resolution and describe the importance 
of this observation for physical random number generation. 


Experiment and Results 



Fig. 1: Experimental configuration. We employ a silicon 
avalanche photodiode, which detects individual photon ar¬ 
rivals. This signal is time-delayed and filtered using an FPGA, 
and the output of the filter drives the modulator which in turn 
varies the light incident on the detector, forming a feedback 
loop. 


Figure 1 shows a schematic of our experimental configu¬ 
ration. Our system has a similar architecture to earlier ex¬ 
periments involving optoelectronic feedback loops, but differs 
in that we use a photon-counting detector, whereas previous 
experiments used an analog photodetector. In either case, 
the signal from the detector is time-delayed and filtered, and 
the output of the filter drives the Mach-Zehnder electro-optic 
modulator (MZM) which in turn controls the light incident on 
the detector, forming a feedback loop. When an analog photo¬ 
diode is used, the feedback loop is modeled by a time-delayed 
nonlinear differential equation. 


-=Ea, + DFm 

I{t) = sin^ — Td) + 


[ 1 ] 


Here, x is the state variable of a linear, time invariant filter, 
matrix E and the vectors F and G describe the characteris¬ 
tics of the filter, I(t) is the normalized intensity of light trans¬ 
mitted through the MZM, and Td is the time delay. When a 
photon-counting detector is used in place of an analog photodi¬ 
ode, the filter variables can be modeled by a linear differential 
equation driven by discrete photon arrivals. In our implemen¬ 
tation, the equations of motion for the filter variables are 


dt 



'T 

0 



where the photon arrivals times, {ti}, are generated by a non¬ 
stationary Poisson point process whose rate, A(t), depends on 
the state of the filter variables. 

A(t) = XoI{t) = Ao sin^ [xi(t - Td) - X 2 {t - Td) + 0] . [3] 
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Fig. 2: Time series, probability distributions, and autocor¬ 
relation functions, a) - c) show experimental data, and d) 
shows the result of a deterministic simulation of Equation (1). 
a) At a low photon rate of XoTd = 12.5, the dynamics ap¬ 
pear Poissonian. The time series has no visible structure, the 
autocorrelation function is sharply peaked at 0, and the dis¬ 
tribution of photon counts in a window of w = Td/4 is nearly 
Poisson, b) XoTd = 200, a slow modulation of the photon rate 
is evident, c) at XoTd = 3200, the photon rate varies smoothly, 
the photon count distribution is bimodal and much wider than 
a Poisson distribution with the same mean, and the autocor¬ 
relation function shows slow oscillations. The deterministic 
simulation d) shows the same features as the high photon rate 
data shown in c). 


In the limit that the Ao is large, the stochastic term in equation 
§ can be replaced with its expectation value, /(t), leading to 
equation 0- 

In our implementation, the time delay is Td=1.734 ms, the 
modulator bias is 0 = 7r/4, and the filter constants Ti = 1.2 
ms, and T2=60 fis. The filter and time delay are implemented 
digitally using an Altera Cyclone II field programmable gate 
array (FPGA) and a digital to analog converter (DAG). The 
clock speed of this device is 151.1515 MHz, and we record 
all of the photon arrival times to this precision. The light 
source in our experiment is a continuous wave fiber-coupled 
distributed feedback laser with a wavelength of 850 nm. Our 
detector has a dark count rate of ~ 100 counts/s and a dead 
time of about 40 ns. We vary the photon rate over a factor 
of 256, from XoTd = 12.5 (7.20 x 10® count/s) to XoTd = 3200 
(1.845 X 10® counts/s). In all of the experiments shown here, 
(3 is kept constant so that fiTd = 8.87. Figure 2 shows several 
time series recorded with this system with increasing photon 
rate, showing a transition from Poisson noise to deterministic 
chaos. We plot Nw{t), the number of photon arrivals in the 
interval [t — re,t]. In Figure 2, all of the plots were generated 
with w = Td/4:. When the incident photon rate is AoTk = 12.5, 
the photons appear to arrive at random, uncorrelated times 
as in a stationary Poisson process. Increasing the incident 
photon rate to XoTd = 200, a smooth modulation of the pho¬ 
ton rate starts to become apparent. At AqT^ = 3200, Nw{t) 
has a smooth character, and qualitatively resembles a low- 
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Fig. 3: Poincare sections. We visualize the emergence of 
a chaotic attractor from Poisson noise with increasing photon 
rate by embedding photon count time series in two dimensions 
with a time delay of A = Td/4, and reducing the dimensional¬ 
ity of the dynamics by plotting points only when the state of 
the system in phase space passes though a codimension-1 sur¬ 
face dehned by xi — X2 = tt. a) - c) show experimental data, 
and d) shows the result of a deterministic simulation. These 
histograms are constructed with a bin width of 1 photon in in 
a) and b), 4 photons in c), and 0.005 in d) 


dimensional chaotic signal. We also p^t the results of a de¬ 
terministic simulation using equation (fy . This time series was 
smoothed with a moving average over a time window of width 
w to be directly comparable with {t)_. We plot the autocor¬ 
relation function, C{t') — — Nw){Nw{t — t') — 

normalized so that the value of the autocorrelation function is 
unity at = 0. As the photon rate increases from XoTd = 12.5, 
the autocorrelation function changes from a (5-like peak, char¬ 
acteristic of a Poisson process, to an oscillatory function which 
shows correlations at long timescales (tens of milliseconds). 
The autocorrelation function of the deterministic simulation 
time series is in close agreement with the autocorrelation func¬ 
tion of the photon arrivals with XoTd = 3200. Histograms of 
Nw (t) also show a transition from a nearly Poisson distribution 
to a bimodal distribution characteristic of the deterministic 
chaotic process. 

To visualize the development of chaos with increasing pho¬ 
ton rate, we show Poincare surfaces of section in Figure 3. 
We perform a time-delay embedding of the experimental time 
series Nw{t), using a time delay of A = Td/4, by construct¬ 
ing a list of points in two dimensional space of the form 
[Nw{t), Nw{t — A)]. Because the attractor has a dimension 
higher than 2, we reduce the dimensionality of the attractor by 
plotting the points only when the state variables pass through 
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Fig. 4: Experimental dependence of growth rate of variance 
on counting time window, w. To indicate the timescale of the 
deterministic dynamics we indicate Td- We also show a line 
indicating Var(iV^u;) = 1, which roughly separates timescales 
over which less than one photon arrives in the counting win¬ 
dow from timescales over which many photons arrive. Distinct 
asymptotic values of Ya,T(Nw)/w are seen in the limits of small 
and large w. The offset between these two values reffects de¬ 
terministic correlations in the photon arrival rate and grows 
with increasing photon rate. 


a codimension 1 Poincare surface dehned by xi—X 2 = tt. The 
embeddings show a similar trend to the plots in Figure 2. We 
see a development of complex chaotic dynamics from discrete 
photon noise as the photon rate increases. The deterministic 
simulation is plotted for comparison, and, as in Figure 2, a 
moving average of width w is employed so that the smoothed 
intensity time series, Iw{t), is directly analogous to Nw{t). 
The deterministic signal in Figure 3d can be regarded as the 
inhnite photon rate limit of the photon counting system. 

Figure 4 shows the dependence of the variance of Nw on 
the window w and offers another indication of the transition 
from shot noise to deterministic chaos. The time integral of an 
uncorrelated random signal executes a random walk in which 
the variance grows linearly with the integration time. For 
this reason, we plot YsiT{Nw)/w in Figure 4. We see distinct 
asymptotic growth rates of the variance with small and large 
w. When w is small, the variance reflects the Poissonian na¬ 
ture of the photon arrivals, and the growth rate of the vari¬ 
ance has roughly constant value of YsiT{Nw)/w = XqL In the 
limit where the counting window is much longer than the time 
scale of the variations in intensity, Nw{t) can be regarded as 
the sum of the photon counts in many independent identically 
distributed intervals, and the central limit theorem implies 
that the variance will grow in proportion to w. As we increase 
the photon rate from XoTd = 12.5 to XoTd = 3200, we see an 
increasing offset between the two asymptotic rates of growth 
of the variance. The variance can be related to the photon 
rate, counting window, and the unnormalized autocorrelation 
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function, ci{t') — ((/(t) — — t') — l)) [^. 


Var(iV^(;) = w 


\qI + Aq 



dt' 



[4] 


Q{w) 


The second term in Equation @ accounts for the difference 
between the observed variance and the variance of a Poisson 
process with the same rate. The quantity 0(ie) has units 
of time, and measures the correlations in I{t) introduced by 
the feedback. This quantity increases from 0 to an asymp¬ 
totic value ©oo as re increases, accounting for the shape of the 
curves shown in Figure 4. In deterministic simulations, we find 
©oo = 150 //s. The value of ©oo is related to the size of the 
intensity fluctuations, and the rate at which G(w) approaches 
this asymptotic value is determined by the timescales of the 
correlations of I{t). 


Entropy quantification. We characterize the entrop y p roduc- 
tion using the {s,t) entropy per unit time, h{e,r) [^. This 
measurement of entropy has two parameters: samjne resolu¬ 
tion, £, and sampling time interval r, which are natural pa¬ 
rameters for most experiments because measurement devices 
record data to finite resolution at discrete times. In addition 
to being experimentally relevant, the dependence of h{e, r) on 
these parameters can reflec t the underlying physical origin of 

unpredictability ME]- . 

In chaotic systems, unpredictability is due to the sensitive 
dependence on initial conditions. Because small perturbations 
grow exponentially in time, chaotic systems generate informa¬ 
tion. The growth of uncertainty is quantified by the Lyapunov 
exponents fii, and in particular the largest exponent, /xi. Pos¬ 
itive Lyapunov exponents and entropy rate both quantify un¬ 
predictability and there is a close relationship between these 
two quantities. One would expect that if a chaotic system 
is sampled infrequently (r/xi 1), successive samples will 
be uncorrelated because of the growth of uncertainty between 
measurements. On the other hand, if the interval between 
successive samples is small (r/xi <C 1), one expects strong 
correlations between adjacent samples and a reduced entropy 
per sample. Experimental and theoretical work using semi¬ 
conductor lasers has shown that these considerations are cru¬ 
cial to physical random number generation using chaotic dy- 
[491. In the limit that r, e: ^ 0, h{e,T) will 


namics 21 22 


approach a hnite valu e, the Kolmogorov-Sinai (or metric) en¬ 
tropy, hks 50 The metric entropy is related to the 
Lyapunov exponents, by 


“ log(2) E 

We calculated the spectrum of Lyapunov exponents from 
Equation (I) |53||54| . There is only one positive Lyapunov 
exponent with a value of /xi/ log(2) = 345 bits/s. The Kaplan- 
Yorke dimension calculated from the Lyapunov spectrum 
is 3.56. 

The (e, r) entropy will have qualitatively different behavior 
as £ ^ 0 depending on the physical origin of unpredictability. 
In chaotic systems, the entropy rate does not depend on either 
the sampling rate or the sampling resolution. This property 
of chaotic systems imposes a theoretical limitation on physical 
random number generation. Increasing the speed and resolu¬ 
tion of a measurement device cannot in principle increase the 
entropy that can be harvested from a deterministic chaotic 


system beyond hks- In contrast to deterministic systems, the, 
entropy rate of stochastic signals diverges like — log( 5 ) for fi¬ 
nite T [38|[40] . 

Another advantage of the {s, r) entropy is that it can be cal¬ 
culated from experiment al da ta using an algorithm described 
by Cohen and Procaccia [^. In our case, we chose to calcu¬ 
late the entropy from NyjJtJ with a counting time window of 
w = Td/4. With this window, Nw{t) approximates the behav¬ 
ior of the deterministic signal as seen in Figures 2 and 3. 


a) 




Fig. 5: {s,t) entropy in deterministic simulation and experi¬ 
ment. Calculations with embedding dimensions of 5,10,15,20 
and 25 are shown, a) Deterministic simulation. The entropy 
is calculated from the intensity, /(t), which is normalized so 
that its values are between I (complete transmission through 
the modulator) and 0. s is measured in the same units as 
I{t). Characteristic of deterministic systems, the entropy rate 
is independent of £, and approaches the largest Lyapunov ex¬ 
ponent as d increases, b) Entropy rate of experimental time 
series. The entropy is calculated from Nw{t), and £ is mea¬ 
sured in photons. At low photon rates we see a divergence 
characteristic of a Poisson processes. As the photon rate in¬ 
creases, the dependence of the entropy on e becomes progres¬ 
sively flatter, and approaches hks- Across photon rates, we 
see a divergence of the entropy rate for small £. The Poisson 
curves were calculated by approximating the Shannon entropy 
of a Poisson process by an integral oyer a Gaussian distribu¬ 
tion with a mean and variance of XqIw. 
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We do not employ an averaging time window to compute the 
entropy from deterministic simulations. 

The first step to computing the entropy rate of an experi¬ 
mental signal is to generate a list of points in d-dimensional 
space using time delay embedding with a delay of r. These 
vectors can be regarded as samples of a d-dimensional proba¬ 
bility distribution over phase space. The entropy of this prob¬ 
ability distribution, Hd, is sometime s re ferred to as the pat¬ 
tern entropy for patterns of length d . In principle Hd can 
be calculated by building a histogram with boxes of width e 
and applying Shannon’s formula, H = —^iPi^og 2 Pi [^. In 
practice, direct application of this approach requires a very 
large amount of data when the embedding dimension is large. 
Cohen and Procaccia propos ed a more efficient algorithm to 
estimate the pattern entropy in the context of estimating 
metric entropy from experimental data. First, one randomly 
selects a small number M of reference points from the time 
series. In our case, M = 5000 was sufficient. For each refer¬ 
ence point z, one computes ni(e), the fraction of points within 
a box of width £ centered on the reference point. The only dif¬ 
ference between a direct calculation of the Shannon entropy, 
and the Cohen-Procaccia procedure is that in a direct calcu¬ 
lation, a rectangular array of bins is used, rather than a set of 
bins centered on random points chosen from the data set. In 
searching for neighbors for the zth reference point, we exclude 
points within a time window of r of that point, as suggested 
by Theiler [^. The pattern entropy is then estimated by 

1 ^ 

Hd{s) = [ 6 ] 

i=l 

It is a general feature of unpredictable signals that Hd grows 
linearly with d in the limit that d is large, and the entropy 
rate is the slope of this linear increase. 

h{e,T) = - lim [Hdis) - Hd-iis)] [ 7 ] 

T d^oo 

Figure 5 shows the entropy per unit time in both determinis¬ 
tic simulation and experiment with r = (3/4)Td. The duration 
of the simulation was 512 x 10® Td- Figure 5a shows that in 
the deterministic simulation, the entropy rate remains flat as 
e decreases. As d increases, this plateau approaches hks, as in¬ 
dicated by Equation In Figure 5b, we see that as the pho¬ 
ton rate increases, the aependence of the entropy on £ becomes 
progressively flatter at high £. Furthermore, in the region that 
this flattening is present, the value of h{£,T) is close to hks- 
The flattening of h{£, r) at high photon rate is another indi¬ 
cation that this system behaves more deterministically in this 
regime. At all photon rates, we see h(£,r) sharply increases 
as e decreases, which is due to the shot noise inherent in the 
system. It is natural to compare the entropy rates we observe 
to a constant-rate Poisson process with the same average rate. 
The Poisson curves in Figure 5 were calculated by approxi¬ 
mating the Shannon entropy of a Poisson distribution by an 
integral over a Gaussian distribution with a mean and vari¬ 
ance of XqIw. In the limit that e ^ \/ XqIw, this leads to the 
asymptotic expression h(£,T) — (—l/r)\og 2 [s/ \/ 27TeXoIw ^, 

indicated by the dashed curves in Fig. 5(b). For small £, 
h{£, r) increases logarithmically with decreasing £, and par¬ 
allels this curve. This logarithmic dependence is more pro¬ 
nounced at lower Aq. 

Discussion 

We show in this paper that the choice of the resolution with 
which we observe our system allows us to see either noisy 
or deterministic dynamics. By counting photon arrivals over 


timescales on the order of the delay time and hlter time 
constants, we see deterministic dynamics in the time series, 
Poincare sections, and the autocorrelation functions. Further¬ 
more, when we observe the dynamics on large scales of both 
value (e) and time {w and r), we find that the entropy rate is 
close in value to the metric entropy calculated from the pos¬ 
itive Lyapunov exponents of the deterministic model, which 
shows that the entropy generation is dominated by the deter¬ 
ministic exponential amplification of small perturbations in 
this regime. 

In contrast, by employing high resolution in photon counts 
and time scales, we see that both the entropy rate and vari¬ 
ance reflect the stochastic nature of the photon arrivals. For 
small values of w, the variance of the number of photon counts 
is equal to the average number of counts, characteristic of a 
Poisson process. The logarithmic dependence of the entropy 
on £ shown in Figure 5 offers another indication of the noisy 
nature of the dynamics at small scales. In addition to showing 
both shot noise and chaos at different scales, our experiment 
also shows a transition from shot noise to chaos with increas¬ 
ing photon rate. The precise control over the rate of photon 
arrivals and dynamical timescales afforded by our experiment 
allows for experimental observation of the interplay of noise 
and dynamics. Our results can be seen to bridge two widely 
used methods of physical random number generation. 

Two prevalent methods have attracted attention for opti¬ 
cal random number generation: those based on si ngle ph oton 
detection from strongly attenuated light sources [59[|60| , and 
those base d on digitized high-speed fluctuations from chaotic 
lasers [^. In the former case, the entropy is claimed to orig¬ 
inate entirely from quantum mechanical uncertainty, yet in 
practice these methods are also subject to unpredictable drift 
and environmental variations. In the latter case, the entropy is 
attributed to the dynamical unpredictability of chaos, but the 
unavoidable presence of spontaneous emission is tho ught to 
play a role in seeding these macroscopic fluctuations [^[^ . 
The system presented here is unprecedented in that it can 
approach macroscopic chaos from the single photon limit, 
thereby revealing the transition from noise to chaos. More¬ 
over, the analysis offers a unified measure of entropy that cap¬ 
tures both behaviors, and clarifies the relationship between 
sampling frequency, measurement resolution and entropy rate. 

The designer of a physical random number generator must 
choose the sampling rate and resolution that they will use 
to collect numbers from a physical system. These decisions 
will impact the entropy rate. Heuristically, finer discretiza¬ 
tion (smaller s) and more frequent sampling (smaller r) lead 
to higher entropy rates, but without the methods presented 
here it is difficult to assess the dependence of the entropy 
rate on these parameters in any given system. The statisti¬ 
cal tests that are usually used to evaluate physical random 
number generation were not designed to answer these 

questions, but rather to certify that a stream of bits is free 
of bias and correlation. If a random number generator em¬ 
ploys post processing (as most do), existing statistical tests 
applied to the output binary sequence provide no insight into 
whether the entropy originates from the physical process or 
the post-processing algorithm employed. The {£,t) entropy 
clarifies the origin and nature of uncertainty and informs the 
choice of sampling rate and measurement resolution. 
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